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Understanding  of  both  water  management  in  PEM  fuel  cells  and  degradation  mechanisms  of  the  gas 
diffusion  layer  (GDL)  and  their  mutual  impact  is  still  at  least  incomplete.  Different  modelling  approaches 
contribute  to  gain  deeper  insight  into  the  processes  occurring  during  fuel  cell  operation.  Considering 
the  GDL,  the  models  can  help  to  obtain  information  about  the  distribution  of  liquid  water  within  the 
material.  Especially,  flooded  regions  can  be  identified,  and  the  water  distribution  can  be  linked  to  the 
system  geometry.  Employed  for  material  development,  this  information  can  help  to  increase  the  life 
time  of  the  GDL  as  a  fuel  cell  component  and  the  fuel  cell  as  the  entire  system.  The  Monte  Carlo  (MC) 
model  presented  here  helps  to  simulate  and  analyse  the  water  household  in  PEM  fuel  cell  GDLs.  This 
model  comprises  a  three-dimensional,  voxel-based  representation  of  the  GDL  substrate,  a  section  of  the 
flowfield  channel  and  the  corresponding  rib.  Information  on  the  water  distribution  within  the  substrate 
part  of  the  GDL  can  be  estimated. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

One  of  the  main  obstacles  on  the  way  of  polymer  electrolyte 
membrane  fuel  cells  (PEMFCs)  towards  commercial  use  is  their 
durability.  Thus  during  the  last  years  a  significant  amount  of 
research  work  has  been  done  to  investigate  PEMFC  degradation 
mechanisms  and  mitigation  strategies.  Up-to-date  overviews  on 
PEMFC  degradation  are  given  in  [  1  -4].  For  a  long  time  the  focus  has 
been  less  on  GDL  ageing  but  different  studies  concerning  the  life¬ 
time  of  fuel  cells  have  rather  investigated  reaction  layer,  membrane 
and  bipolar  plate  ageing.  Nevertheless,  recent  results  show  that 
GDL  degradation  should  not  be  neglected  anymore.  For  the  case  of 
the  electrode  in  [5]  it  was  demonstrated  that  the  degradation  of  the 
PTFE  results  in  approximately  two  times  higher  performance  losses 
than  the  catalyst  degradation  by  Pt  agglomeration  after  1000  oper¬ 
ating  hours  under  the  operating  conditions  chosen;  this  at  least 
implies  that  an  impaired  water  household  caused  by  the  loss  of 
GDL  hydrophobicity  may  have  similarly  severe  effects.  Though  the 
correlation  between  liquid  water  and  GDL  degradation  has  been 
discussed  since  quite  some  time  [6],  the  mechanisms  are  still  not 
well  understood.  In  principle  two  reasons  of  degradation  for  the 
PEMFC  GDL  are  presumed: 
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•  Changes  of  the  wetting  behaviour  of  the  GDL:  The  wetting 
behaviour  changes  may  be  induced  by  PTFE  decomposition,  see 
e.g.  [7],  or  depletion  caused  by  oxidation  of  the  PTFE  or  the  car¬ 
bon,  mechanical  effects  [8]  (compression  or  fibre  cracking)  or 
impurities  [6].  By  these  effects  the  originally  hydrophobic  GDL 
becomes  more  hydrophilic  which  results  in  an  accumulation  of 
liquid  water  in  the  GDL  pores  and  blocking  of  these  pores  for 
reactant  transport. 

•  Changes  of  the  structure  of  the  GDL  caused  by  mechanical  effects 
(compression  below  the  ribs,  see  [8,9])  or  carbon  oxidation  (for 
the  case  of  the  MEA,  see  e.g.  [10])  could  change  the  electric  as 
well  as  the  heat  conductivity  of  the  GDL. 

These  mechanisms  are  not  limited  just  to  the  GDL,  but  also  influ¬ 
ence  other  components  of  the  PEM  fuel  cell.  This  phenomenon 
demonstrates  how  important  the  interaction  of  the  different  degra¬ 
dation  mechanisms  is.  The  changes  in  the  wetting  behaviour  and 
the  structure  of  the  GDL  and  thus  its  durability  depend  on  the  oper¬ 
ating  conditions  of  the  PEM  fuel  cell.  In  [11]  it  is  stated  that  the 
presence  of  liquid  water,  excess  reactants  and  high  current  densi¬ 
ties  induces  significant  degradation  of  the  MEA  and  the  GDL  The 
authors  of  [11,12]  observe  that  high  humidification  and  especially 
liquid  water  are  correlated  with  PTFE  decomposition. 

The  task  of  modelling  the  GDL  water  household  implies  two 
challenges:  finding  an  appropriate  description  of  the  porous  struc¬ 
ture  and  modelling  the  dynamic  behaviour  of  water  and  gas  within 
this  structure.  In  literature,  different  modelling  approaches  can  be 
found  associated  with  porous  media  description  and  modelling  of 
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transport  phenomena  in  porous  media.  On  the  one  hand  there  are 
several  approaches  to  find  a  reasonable  description  of  the  porous 
media,  on  the  other  hand  there  are  approaches  to  model  the  trans¬ 
port  processes  in  porous  media  by  solving  macroscopic  transport 
equations  or  by  describing  molecular  particle  transport.  The  for¬ 
mer  can  be  divided  into  models  representing  the  real  geometry  of 
the  porous  structure,  statistically  generated  structures,  randomly 
generated  structures,  abstracted  structures  representing  certain 
properties  of  the  real  material;  even  continuous  material  linked 
with  some  integral  structural  parameters  as  permeability  can  be 
used  to  describe  transport  properties  of  porous  structures. 

1.1.  Macroscopic  transport  models 

As  an  example  for  a  macroscopic  modelling  approach,  compu¬ 
tational  fluid  dynamics  (CFD)  can  be  mentioned.  This  approach 
can  be  subdivided  into  two  major  categories.  The  first  category 
is  considering  the  real  multi-fluid  approach  where  the  entire  set 
of  conservation  equations  for  each  phase  is  solved  and  the  cou¬ 
pling  of  the  equations  is  achieved  through  pressure  and  inter-phase 
exchange  coefficients.  A  successful  example  of  an  implementation 
can  be  found  in  [13].  The  second  category  is  the  multi-phase  mix¬ 
ture  approach  where  the  mixture  momentum  equations  are  solved 
and  the  relative  velocities  are  calculated  in  a  post-iterative  step. 
One  of  the  first  papers  introducing  this  approach  is  [14].  In  the 
context  of  CFD,  the  porous  material  can  be  described  both  by  real 
geometry  data  or  continuously  together  with  integral  structural 
parameters. 

A  further  approach  which  has  been  successfully  applied  to  the 
field  is  the  Lattice  Boltzmann  (LB)  method  which  indirectly  solves 
the  Navier-Stokes  equation  by  dealing  with  the  Boltzmann  equa¬ 
tion  for  fluid  particle  movement  [  1 5  ].  It  works  on  a  lattice  consisting 
of  equally  shaped  unit  volumes  (so-called  lattice  sites,  or  cells), 
which  are  either  solid  or  filled  with  fluid. 

Pore  network  models  are  designed  to  capture  the  transport 
properties,  but  not  the  morphology  of  porous  media,  i.e.  they  use 
abstracted  structures  to  represent  the  properties  of  the  pore  space. 
The  use  of  networks  of  capillary  tubes  to  model  porous  media  was 
first  suggested  in  [16]. 

1.2.  Nano-  and  microscopic  models 

Models  based  on  macroscopic  liquid  properties  may  be  use¬ 
ful  up  to  at  least  the  p,m  scale.  Going  to  even  smaller  structures, 
e.g.  the  nm-scale  pores  within  the  MPL  material  of  a  fuel  cell 
GDL,  the  model  assumptions  break  down  and  fluid  flow  must  be 
described  on  the  basis  of  appropriate  particle  interactions.  Molec¬ 
ular  dynamics  (MD)  modelling  is  one  of  the  principal  tools  to 
theoretically  study  fluid  flow  based  on  molecular  interactions.  This 
computational  method  describes  the  time  dependent  behaviour  of 
a  molecular  system,  thereby  calculating  particle  interactions  either 
using  empirical  forcefields  or  even  quantum  mechanical  methods. 
Commonly,  to  describe  particle  movement,  the  classical  equations 
of  motion  are  integrated  based  on  an  adequate  discretization  in 
time.  The  MD  method  has  been  routinely  used  to  investigate  the 
structure,  dynamics  and  thermodynamics  of  fluids.  With  MD  sim¬ 
ulations  one  can,  e.g.  study  both  thermodynamic  properties  and/or 
time  dependent  (kinetic)  phenomena  such  as  fluid  flow  in  porous 
material  (for  a  review,  see  [17-20]). 

A  further  method  which  has  been  commonly  used  on  the 
molecular  scale  is  the  Monte  Carlo  (MC)  method,  which  was  first 
described  in  [21  ].  Since  then,  MC  simulations  have  been  used  in  sev¬ 
eral  fields  of  scientific  modelling.  Concerning  fluid  behaviour,  some 
authors  focus  on  phase  change  as  [22  ]  or  on  two  phase  flow  without 
phase  change  as  for  example  [23].  But  also  percolation  theory  [24] 
and  diffusion,  including  Knudsen  diffusion,  [25]  are  considered.  For 


all  of  these  models,  a  proper  representation  mirroring  the  structure 
of  the  porous  medium  has  to  be  found  first.  Therefore,  depending 
upon  the  physical  model  to  be  used,  simple  2D  cylindrical  pores  as 
in  [22]  or  spherical  clusters  as  in  [25]  are  employed. 

Furthermore,  the  MC  method  has  been  used  to  gain  a  represen¬ 
tation  of  the  structure  of  porous  media  and  the  associated  effective 
properties  (see  [22-29]).  Some  of  the  MC  models  also  focus  on 
the  determination  of  general  parameters  of  porous  media,  e.g.  in 
[28]  and  [29],  the  MC  simulation  is  used  to  compute  the  effective 
diffusivities  in  randomly  overlapping  fibres  and  filament  bundles 
representing  a  porous  medium.  Further  articles  especially  focus  on 
fuel  cell  topics:  in  [26],  electrical  conduction  aspects  are  examined, 
in  [27],  the  MC  method  is  used  to  generate  the  geometry  for  a  lattice 
Boltzmann  calculation  of  the  permeability. 

There  are  further  particle-based  models  using  not  explicitly 
molecular,  but  small-scale  descriptions  of  material  behaviour. 
These  include  randomly  distributed,  overlapping  grains  such  as 
oblate  and  prolate  spheroids.  They  are  known  to  quite  accurately 
describe  the  pore  space  in  sandstone,  foams  and  fibre  networks 
[30-33]. 

The  novel  approach  presented  in  this  article  employs  the  MC 
method  not  on  the  molecular,  but  on  the  p,m  scale  to  describe  water 
behaviour  within  randomly  generated  porous  material  on  a  statis¬ 
tical  basis,  i.e.  without  the  need  of  numerically  solving  transport 
equations.  Also  real  geometry  data  could  be  used. 

1.3.  Experimental  visualization  of  water  in  PEMFC  GDLs 

Experimental  in  situ  visualization  of  liquid  water  within  PEM- 
FCs  has  been  approached  using,  e.g.  optical  imaging  [34],  neutron 
radiography  and  tomography  [35],  synchrotron  radiography  or 
tomography.  Combined  optical  imaging  and  neutron  radiography 
studies  have  been  performed  [36],  but  the  spatial  resolution  of 
this  technique  is  too  low  to  gain  insight  into  inner-GDL  processes. 
Employing  synchrotron  radiography,  e.g.  Manke  et  al.  [37]  and 
Hartnig  et  al.  [38]  have  shown  that  in  situ  water  distributions  within 
fuel  cell  GDLs  of  appropriate  resolution  can  also  be  obtained  in 
through-plane  direction.  Unfortunately,  to  our  knowledge,  so  far 
no  according  experimental  data  dealing  with  ageing  or  variation  of 
PTFE  coverage  is  available. 

Clearly,  the  experimental  progress  on  the  field  opens  new 
perspectives  concerning  water  household  analysis.  Nevertheless, 
further  effort  on  the  field  of  theoretical  modelling  is  required  to 
provide  and  facilitate  deeper  understanding  and  to  relieve  material 
improvement  by  predicting  the  consequences  of  changing  param¬ 
eters.  The  MC  model  described  in  the  following  is  intended  to 
complement  the  established  tools  in  this  field. 

This  article  is  structured  in  the  following  way:  first  the  MC 
model  developed  to  study  the  water  distribution  in  fuel  cells  will  be 
described.  Also  model  specific  properties  as  the  handling  of  porous 
media  and  the  steady  state  as  approximation  of  real  fuel  cell  con¬ 
ditions  will  be  presented.  The  results  part  is  discussing  fuel  cell 
systems  of  different  PTFE  content  which  represent  different  ageing 
situations.  Several  approaches  to  analyse  the  simulations  are  pre¬ 
sented  and  discussed.  The  article  will  be  completed  by  a  conclusions 
and  outlook  section. 

2.  Monte  Carlo  model 

2.1.  Description  of  the  MC  model 

In  contrast  to  the  commonly  known  MC  models,  the  model  pre¬ 
sented  in  the  following  does  not  work  on  the  molecular  level,  but 
on  the  fxm  scale.  MC  models  working  on  the  molecular  scale  often 
use  a  (N,  V,  T)-ensemble  to  simulate  a  system  at  equilibrium.  This 
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means  that  using  a  constant  number  of  particles  N,  volume  V  and 
temperature  T,  the  system  explores  the  microscopic  phase  space  in 
a  way  that  average  values  gained  from  the  simulation  correspond  to 
the  values  for  the  macroscopic  thermodynamic  equilibrium  state. 

The  model  presented  here  is  grid-based  with  a  voxel  size  of 
5  [xm  x  5  [xm  x  5  [xm.  The  lattice  sites  can  be  allocated  with 

•  fibre  material:  graphite  or  PTFE; 

•  liquid  water; 

•  “nothing”  (gas  phase  not  explicitly  included); 

Within  the  model  periodic  boundary  conditions  are  applied  in 
all  directions  of  space.  These  boundary  conditions  can  be  lifted  by 
walls  delimiting  the  system.  The  wall  particles  do  not  have  any 
interaction  to  the  water  droplets. 

During  a  MC  iteration,  each  water-occupied  voxel  is  chosen  once 
in  random  order  for  a  tentative  move,  referred  to  as  single  MC  step 
in  the  following.  One  free  position  among  the  six  next  neighbour 
voxels  is  randomly  selected  and  the  energy  of  the  new  and  the  old 
state  are  calculated.  These  interactions  between  the  occupied  vox¬ 
els  are  modelled  as  pair  interactions.  The  values  of  the  interactions 
are  based  on  surface  energies.  For  the  water-water  interactions, 
-72mNnrr1  at  20°C  are  used.  Temperature  dependence  of  the 
value  is  taken  into  account  via 

Eh2o  =  0.07275  (l  -  0.002  (^  -  29l))  . 

For  graphite-water  -95  mN  irr 1  and  for  PTFE-water  +30  mN  m-1 
have  been  estimated.  The  energy  of  a  given  state  is  calculated  tak¬ 
ing  into  account  all  occupied  neighbouring  positions  of  the  voxel 
under  consideration.  Direct  neighbours  are  fully  counted,  whereas 
for  diagonal  neighbours  a  weighting  factor  of  1/r6  is  applied.  The 
resulting  energy  is  converted  into  a  molar  energy. 

As  usually,  the  total  energy  difference  AF  of  the  possible  new 
energy  to  the  old  state  determines  the  acceptance  probability  for 
the  MC  step  via  the  Boltzmann  factor 


Hereby,  T  is  the  absolute  temperature  and  R  denotes  the  universal 
gas  constant.  If  the  energy  difference  AF  has  a  negative  value,  the 
factor  is  bigger  than  1  and  the  MC  move  is  accepted.  If  AF  is  posi¬ 
tive,  the  exponent  is  smaller  than  1  and  a  random  number  between 
0  and  1  is  calculated.  If  the  Boltzmann  factor  is  larger  than  this  ran¬ 
dom  number,  the  MC  move  is  accepted.  Otherwise,  the  old  state 
of  the  system  is  preserved.  In  both  cases,  the  iteration  continues 
with  the  random  selection  of  the  next  water-occupied  voxel.  The 
pseudo  random  number  generator  used  within  our  model  is  the 
MT19937  ‘Mersenne  Twister’  from  the  Gnu  Scientific  Library  (GSL) 
(for  further  information,  see  [39]). 

2.2.  Steady  state 

An  operating  fuel  cell  cannot  be  adequately  described  by  a 
system  in  thermodynamic  equilibrium.  Hence,  real  FC  operating 
conditions  need  to  be  approximated  by  other  means.  FC  operation 
at  constant  current  and  flow  rates  can  be  characterized  as  a  steady 
state.  But  as  within  the  framework  of  a  MC  simulation,  there  is 
no  meaningfull  definition  of  the  term  ‘time’,  also  the  concepts  of 
current  and  flow  cannot  be  directly  applied.  Therefore,  we  have 
developed  a  method  to  mimic  steady  state  conditions  within  a  MC 
simulation. 

Real  cell  operating  conditions  employing  constant  current  and 
flow  rates  correspond  to  a  constant  amount  of  water  being  trans¬ 
ported  through  the  GDL  and  then  being  removed  by  the  gas  flow 
through  the  channel  structure  within  a  given  time  span.  Within  the 


context  of  our  model,  this  situation  can  be  represented  by  enforc¬ 
ing  a  certain  number  of  water  ‘particles’  to  be  created  (source)  or 
destroyed  (sink),  respectively,  in  different  regions  of  the  system. 
Hereby,  the  aim  is  to  reach  a  state  in  which  the  average  number  of 
water  particles  created  and  destroyed  during  a  MC  iteration  reaches 
an  equal  and  constant  value.  To  describe  a  certain  number  of  events 
per  iteration ,  the  term  ‘rate’  is  used  in  the  following,  although  we 
are  aware  that  this  implies  a  re-definition  of  an  established  time- 
related  concept  which  bears  the  danger  of  causing  confusion. 

To  simulate  a  cathodic  GDL  within  a  FC  in  operation,  the  water 
source  is  located  within  the  MPL  to  mimic  water  production  at  the 
reaction  layer.  Water  creation  is  again  a  mainly  random  process 
taking  place  at  free  voxels  within  the  MPL  pores,  governed  by  a 
certain  creation  probability.  Nevertheless  it  is  influenced  to  some 
extend  by  the  water  coverage  of  the  surroundings.  Water  annihi¬ 
lation  only  takes  place  within  the  channel  to  mimic  water  removal 
by  the  gas  flow  there.  Only  water  particles  having  a  surface  towards 
free  space  can  be  annihilated,  again  governed  by  a  certain  predeter¬ 
mined  probability.  Both  predetermined  probability  values  do  not 
have  a  direct  meaning,  many  combinations  of  values  may  lead  to 
the  same  converged  rate. 

Besides  the  requirement  of  reaching  such  a  converged  source 
and  sink  rate,  the  usual  condition  for  MC  simulations  that  the  mean 
number  of  accepted  MC  steps  per  iteration  should  be  converged  is 
taken  as  a  further  criterion.  Additionally,  the  fraction  of  pore  space 
in  the  inner-GDL  region  which  is  occupied  by  water  is  required  to 
be  constant.  This  quantity  will  be  denoted  as  inner-GDL  hydration 
in  the  following.  If  all  three  conditions  are  fulfilled,  the  state  of  the 
system  shall  be  denoted  as  steady-state. 

2.3.  Handling  of  porous  structures 

Fig.  1  illustrates  the  part  of  the  fuel  cell  explicitly  represented 
within  the  MC  model.  As  can  be  seen  from  the  picture,  not  the  com¬ 
plete  fuel  cell  is  modelled,  but  only  the  cathode  side  including  the 
channel,  the  GDL  and  parts  of  the  ribs  of  a  flowfield. 

The  GDL  fibre  structure  can  be  generated  within  the  model 
employing  an  algorithm  based  on  random  growth  processes.  GDL 
fibres  with  a  thickness  of  10  |xm  are  grown  in  random  directions 
until  a  maximum  length  or  an  obstacle  is  reached.  Consequently, 
fibre  intersections  are  avoided.  Fibres  are  generated  until  a  given 
porosity  value  is  reached  (porosity  value:  76%  within  this  study). 
For  the  GDL  calculations  presented  in  the  following,  a  ‘paper-like’ 
structure  is  generated  by  preferring  fibre  growth  within  the  xy- 
plane  over  the  z-direction.  Additionally,  a  MPL  model  consisting  of 
hydrophobic  (PTFE)  material  of  arbitrary  density  is  added. 

As  an  additional  feature,  an  interface  to  enable  both  importing 
of  an  externally  generated  structure,  e.g.  originating  from  ade¬ 
quate  postprocessing  of  tomography  data  (e.g.  reconstructed  from 
synchrotron  tomography  data),  and  reading  of  program  generated 
structures  has  been  defined.  The  latter  feature  can  be  used  to  sys- 
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tematise  simulation  studies  by  using  the  same  basic  structure  for 
all  runs.  This  enables  to  compare  model  runs  with  different  PTFE 
content  as  described  in  the  following  Section  3. 

As  the  simulation  box  only  covers  a  small  section  within  the 
flowfield,  periodic  boundary  conditions  are  applied  in  x  andy  direc¬ 
tions,  whereas  the  system  is  delimited  by  a  wall  in  z-direction. 
Though,  simulating  realistic  dimensions  for  the  channel  would 
impose  too  high  computational  effort,  so  the  dimensions  in  x  and 
y  directions  have  to  be  reduced,  whereas  thickness  has  a  real¬ 
istic  value.  The  following  case  study  always  uses  a  system  size 
of  50  x  50  x  58  voxels,  i.e.  250  p,m  x  250  p,m  x  290  p,m.  Material 
compression  below  the  ribs  can  also  be  reflected,  a  value  of  20% 
compression  is  applied  in  the  following. 

3.  Results 

3.1.  Simple  systems  in  equilibrium 

For  testing  and  validating  the  MC  model,  a  few  simple  test  cases 
were  studied  at  constant  N,  V  and  T.  At  first,  a  random  water  dis¬ 
tribution  in  free  space  with  zero  gravity  was  generated.  Gravity 
can  be  neglected  due  to  the  fact  that  the  surface  energy  values 
are  by  several  orders  of  magnitude  larger  than  gravity.  When  the 
equilibrium  state  is  reached,  a  single  water  drop  is  formed;  this 
geometry  minimizes  the  surface  energy  and  is  the  one  which  should 
be  expected  from  the  physics  of  the  system.  Furthermore  water 
behaviour  on  hydrophilic  and  hydrophobic  material  was  studied. 
In  equilibrium  the  water  drop  adsorbs  to  the  preferable,  hydrophilic 
material  and  forms  a  contact  angle  below  90°  to  the  solid  sur¬ 
face.  In  summary,  the  basic  physical  behaviour  of  liquid  water  as 
drop  formation  and  behaviour  against  hydrophilic  and  hydrophobic 
material  is  reflected  correctly. 

3.2.  GDI  ageing 

As  mentioned  in  Section  1,  two  mechanisms  are  deemed  to  be 
responsible  for  GDL  degradation:  changes  of  the  wetting  behaviour 
due  to  PTFE  loss  and  structural  changes  caused  by  mechanical 
effects  (e.g.  compression  below  the  ribs).  To  model  the  effect  of 
PTFE  loss,  several  systems  with  decreasing  PTFE  content  were  set 
up,  namely  85%  PTFE  coverage  as  model  for  a  new  GDL,  and  75%, 
65%  and  55%  PTFE  coverage  to  represent  different  stages  of  ageing. 
PTFE  coverage  in  our  case  means  the  fraction  of  fibre  surface  cov¬ 
ered  by  PTFE,  not  the  fraction  of  the  total  weight.  The  latter  value  is 
sometimes  given  by  GDL  manufacturers,  but  the  direct  relation  to 
the  surface  coverage  remains  unclear.  Thus  we  have  chosen  to  take 
fraction  of  PTFE  covered  surface  as  input  parameter.  As  a  general 
scheme  to  set  up  a  GDL  structure  for  different  MC  simulation  runs, 
first  the  pure  carbon  ‘backbone’  of  fibres  is  generated  (see  Fig.  2(a)). 
The  resulting  setup  constitutes  the  basic  structure  which  can  be  re¬ 
used  for  all  simulations.  Fig.  2(b)  and  (c)  show  a  85%  PTFE  coverage 
and  a  55%  PTFE  coverage  for  the  given  backbone  structure.  For  all 
runs  within  the  study,  the  number  of  MC  iterations  after  reach¬ 
ing  the  steady  state,  i.e.  those  which  are  used  for  analysis,  is  about 
3  x  106.  The  converged  source  and  sink  rates  are  reaching  the  same 
value  in  all  cases. 

Snapshots  of  the  resulting  water  distribution  are  shown  in  Fig. 
3(a)  and  (b)  for  PTFE  contents  of  85%  and  55%.  In  the  case  of  85% 
PTFE  content,  one  can  only  find  a  few  water  drops  (represented  by 
blue  color)  within  the  GDL.  Some  water  can  also  be  found  on  the 
channel  walls  and  within  the  channel.  In  contrast,  regarding  the 
snapshot  of  the  GDL  featuring  a  PTFE  coverage  of  55%,  all  free  pores 
of  the  GDL  seem  to  be  filled  with  water  and  in  addition  the  number 
of  water  particles  on  the  channel  walls  and  within  the  channel  has 
increased.  Water  drops  levitating  within  the  channel  should  not 


be  interpreted  as  such.  Although  they  can  be  observed  in  the  sim¬ 
ulation  snapshots,  they  correspond  to  energetically  unfavourable, 
though  possible,  configurations.  LJnder  real  two-phase  flow  con¬ 
ditions,  water  would  be  immediately  removed  by  the  gas  flow 
within  the  channel.  As  can  be  seen  in  Section  3.3,  they  completely 
disappear  during  statistical  averaging,  as  should  be  expected.  Con¬ 
sequently,  they  are  irrelevant  for  the  water  distribution  resulting 
from  the  simulations. 

As  a  first  deeper  analysis,  the  mean  trial  move  acceptance  and 
the  inner  hydration  for  different  PTFE  contents  are  considered.  The 
trial  move  acceptance  is  formed  as  a  mean  value  of  the  accepted  MC 
moves  per  MC  iteration  (see  Fig.  4(a)).  The  higher  the  PTFE  content, 
the  more  moves  are  accepted.  The  number  of  accepted  moves  on 
the  one  hand  depends  on  the  number  of  water  particles  within  the 
GDL,  on  the  other  hand  on  the  material  properties.  The  more  water 
particles  are  present  within  the  GDL,  the  less  MC  moves  are  possible 
because  of  the  lower  probability  to  find  a  free  space  to  move  to.  A 
higher  value  of  the  acceptance  rate  also  means  that  water  particles 
move  more  frequently,  cf.  Section  3.3. 

Fig.  4(b)  represents  the  percentage  of  lattice  sites  occupied  with 
water  within  the  GDL.  As  could  be  expected,  the  water  content 
within  the  inner  GDL  region  strongly  increases  with  decreasing 
PTFE  content.  Notably,  going  from  85%  to  65%,  the  value  roughly 
doubles  for  each  step.  In  contrast,  between  65%  and  55%  teflon  cov¬ 
erage,  an  increase  of  over  300%  can  be  observed.  This  feature  will 
be  analysed  in  detail  in  Section  3.3. 


3.3.  Analysis  of  water  clustering  and  connectivity 

To  gain  further  insight  into  the  structure  of  the  water  phase 
within  the  simulated  GLD  material  and  to  be  able  to  assess  its  pos¬ 
sible  impact  on  fuel  cell  operation,  both  finding  a  measure  for  the 
mean  properties  (i.e.  ensemble  average)  of  the  water  phase  and 
characterizing  its  geometric  structure  are  required.  For  the  first 
task,  one  rather  obvious  approach  is  to  calculate  for  each  voxel  the 
probability  of  being  occupied  by  water  during  the  steady  state  sim¬ 
ulations.  By  defining  a  limit  above  which  a  voxel  should  be  regarded 
as  on  average  occupied  by  water,  a  mean  water  distribution  can  be 
derived. 

Yet,  though  this  enables  visualization  which  helps  to  get  an 
impression  on  how  the  GDL  material  is  being  filled  by  water,  fur¬ 
ther  consideration  of  the  internal  structure  of  the  water  phase  is 
necessary  to  be  able  to  quantify  and  rank  its  properties.  One  inter¬ 
esting  question  is  whether  the  water  is  distributed  approximately 
equally  over  the  available  empty  space,  or  if  it  rather  forms  small  or 
even  big  clusters  and  how  these  are  interconnected.  This  informa¬ 
tion  can  be  used  to  understand  why  routes  of  gas  transport  may  be 
blocked  or,  on  the  other  hand,  water  connections  bridging  the  fibre 
structure  and  forming  routes  of  water  transport  to  the  channel  may 
be  formed. 

The  so-called  Hoshen-Kopelman-algorithm  [40]  provides  a 
commonly  used  tool  to  analyse  cluster  structures.  Its  rationale  is  to 
analyse  whether  particles,  i.e.  water-occupied  voxels  in  this  case, 
have  neighbours  of  the  same  material,  i.e.  likewise  water-occupied 
voxels.  Particles  which  are  connected  in  this  way  are  assigned  to 
belong  to  the  same  cluster,  which  in  the  end  results  in  a  cluster  map 
containing  information  on  the  total  number  of  clusters  in  the  sys¬ 
tem,  the  respective  size  of  each  cluster  and  which  water-occupied 
voxels  belong  to  it. 

To  relate  the  water  clustering  properties  to  the  PTFE  load,  i.e.  the 
presupposed  ageing  effect,  a  Hoshen-Kopelman-analysis  of  several 
snapshots  from  each  steady  state  simulation  for  the  different  PTFE 
coverage  values  has  been  performed.  A  typical  example  for  each 
coverage  is  shown  in  Fig.  5.  To  visualize  the  FIoshen-Kopelman 
results,  each  distinguishable  cluster  is  assigned  a  different  color. 
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(a) 

Fig.  2.  Pure  carbon  ‘backbone’  structure  (a)  and  two  different  PTFE  loadings  (b)  85%  PTFE,  and  (c)  55%  PTFE.  Grey  represents  the  graphite  fibres  and  the  flowfield  ribs  consisting 
of  graphite.  PTFE  covered  voxels  are  displayed  in  green  color.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of 
the  article.) 


(a)  (b) 


Fig.  3.  Water  distribution  for  different  PTFE  loadings.  Blue  represents  the  water,  grey  the  graphite  fibres  and  the  flowfield  ribs  and  green  the  PTFE.  (a)  85%  PTFE  coverage; 
(b)  55%  PTFE  coverage.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


As  can  be  seen  from  Fig.  5(a),  which  represents  a  new  GDL,  the 
overall  quite  low  water  content  within  the  GDL,  cf.  Fig.  4(b),  is 
apparently  randomly  distributed  over  the  pore  space  in  small  clus¬ 
ters.  The  only  big  clusters  that  can  be  identified  constitute  a  water 
film  on  the  hydrophilic  channel  walls.  The  latter  to  some  extend 
reach  under  the  channel  ribs,  where  first  larger  agglomerations  of 
water  clusters  can  be  observed. 

Going  to  the  teflon  coverages  below  85%  modelling  aged  mate¬ 
rial,  as  can  be  expected  one  finds  an  increasing  water  occupation 
with  decreasing  PTFE  content.  At  75%,  especially  the  water  clusters 
residing  under  the  flowfield  ribs  grow  remarkably.  This  trend  is 
affirmed  for  65%  together  with  the  emergence  of  larger  clusters  in 


(a) 


the  inner  regions  of  the  GDL.  Still,  they  remain  separated  and  seem 
randomly  distributed. 

In  contrast  to  the  continuous  increase  so  far,  going  to  55%  brings 
along  an  escalating  water  content,  Fig.  4(b).  Whereas  for  the  other 
ageing  grades  the  inner-GDL  water  content  roughly  doubled  per 
1 0%  points  decrease  of  PTFE,  we  now  see  an  increase  to  above  300%. 
As  Fig.  5(d)  shows  this  already  corresponds  to  a  flooding  of  large 
parts  of  the  inner-GDL  volume  resulting  in  a  water  phase  which 
is  mostly  continuous,  i.e.  mainly  consists  of  one  large  cluster.  Of 
course,  although  this  liquid  water  content  still  should  leave  free 
pore  space  for  the  gas  flow,  such  a  large  cluster  extending  over  the 
whole  GDL  area  already  constitutes  a  major  obstacle.  Furthermore, 


Fig.  4.  The  trial  move  acceptance  over  PTFE  content  (a)  represents  the  percentage  of  accepted  moves.  The  inner  hydration  over  PTFE  content  (b)  indicates  the  percentage  of 
lattice  sites  occupied  with  water. 
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Fig.  5.  Formation  of  water  clusters  for  different  PTFE  content  (multi-colored),  ribs  and  fibres  are  displayed  as  grey  scaffold;  snapshots  of  single  MC  configurations,  (a)  85% 
PTFE;  (b)  75%  PTFE;  (c)  65%  PTFE;  (d)  55%  PTFE. 


as  liquid  water  within  the  GDL  in  turn  is  suspected  to  enhance  GDL 
degradation  [11],  one  may  expect  a  self-enhancing  of  ageing  from 
here  on. 

Of  course,  the  argumentation  presented  so  far  is  based  on  simu¬ 
lation  snapshots,  and  meaningful  conclusions  from  MC  simulations 
should  be  based  on  average  values  representing  the  statistics  of  the 
ensemble.  But  for  the  case  of  water  distribution,  the  latter  are  not 
as  easily  interpretable  as,  for  e.g.  values  of  energy,  etc.  which  are 
commonly  averaged.  It  is  rather  easy  to  define  a  mean  water  distri¬ 
bution,  e.g.  by  defining  that  all  voxels  which  feature  a  occupation 
probability  of  70%  or  higher  over  the  respective  steady-state  sim¬ 
ulation  should  be  regarded  as  water-occupied.  But  on  the  other 
hand,  one  must  bear  in  mind  that  a  water  drop  which  is  alternately 
changing  its  geometric  location  between  two  voxels  would  not  con¬ 
tribute  to  the  distribution  defined  this  way,  although  it  is  present  at 
nearly  the  same  place  in  all  configurations.  On  the  other  hand,  when 
finding  occupation  probabilities  of  over  70%  over  ample  regions  of 
space,  one  can  be  rather  sure,  that  this  corresponds  to  a  significant 
water  agglomeration  that  can  be  located  at  this  spot.  If  furthermore 
the  findings  from  such  mean  values  are  consistent  with  represen¬ 
tative  screenshots  as  presented  before  and  integrated  data  as  the 
distributions  presented  later,  they  can  be  deemed  to  be  meaningful. 

Fig.  6  features  such  water  cluster  data  for  a  occupation  prob¬ 
ability  of  over  70%  as  described  just  above.  Additionally,  again 
a  Hoshen-Kopelman  analysis  has  been  performed  according  to 
which  different  clusters  are  drawn  using  different  colors  as  for 
the  snapshot  cases  depicted  above.  As  mentioned,  by  definition 
this  average  water  occupation  tends  to  contain  only  that  part  of 
the  overall  water  content  of  the  steady-state  simulations  which 


remains  stationary.  This  can  be  clearly  seen  from  a  comparison 
of  Figs.  5(a)  and  6(a):  Though  the  total  water  content  is  constant 
throughout  the  simulation,  the  numerous  small  water  clusters  vis¬ 
ible  in  the  snapshot  picture  have  no  counterpart  in  Fig.  6(a).  Flere, 
one  just  finds  the  water  film  on  the  channels  walls  as  known  from 
above  and  some  water  occupation  below  the  flowfield  ribs.  This 
means  that  most  of  the  liquid  water  content  within  the  GDL  is 
not  only  dispersed  over  the  structure,  but  also  the  resulting  small 
clusters  are  rather  mobile  within  the  simulation,  i.e.  many  config¬ 
urations  of  similar  energy  can  be  realized  by  the  system. 

The  small  pillars  at  the  bottom  part  of  Fig.  6(a),  which  are  also 
visible  in  the  other  figures,  correspond  to  pores  within  the  MPL, 
which  are  used  to  model  water  transport  from  the  reaction  layer 
towards  the  fibre  structure;  this  is  rather  a  technical  feature,  as  due 
to  the  different  pore  size  ranges  of  MPL  and  fibre  structure,  simu¬ 
lating  the  MPL  geometry  on  the  same  lengthscale  is  not  feasible. 
Consequently,  they  should  not  be  over-interpreted. 

Going  towards  75%  and  65%  PTFE  coverage,  obviously  more  of 
the  water  clusters  become  resident;  stationarity  appears  to  be  a 
feature  which  is  related  to  increasing  cluster  size.  Furthermore,  the 
water  agglomerations  below  the  flowfield  ribs  as  monitored  from 
the  snapshot  examples  can  be  confirmed.  Regarding  55%  teflon  cov¬ 
erage,  the  formation  of  one  single  and  -  compared  to  the  other  cases 
-  extraordinary  big  connected  water  cluster  can  be  validated.  Even 
for  the  relatively  demanding  criterion  of  at  least  70%  probability,  it 
continuously  connects  the  MPL  region  to  the  channel  ribs.  Overall, 
one  can  state  that  the  results  of  the  occupation  probability  analysis 
back  and  clarify  the  findings  from  the  snapshot  examples.  Further¬ 
more,  they  enable  the  interpretation  in  terms  of  cluster  mobility. 
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Fig.  6.  Mean  hydration  for  different  PTFE  content  values.  Voxels  featuring  at  least  70%  probability  of  water  occupation  are  shown  as  water,  cluster  affiliation  represented  by 
multi-coloring,  (a)  85%  PTFE;  (b)  75%  PTFE;  (c)  65%  PTFE;  (d)  55%  PTFE. 


This  water  cluster  formation  within  the  MC  model  appears  to  be 
in  accordance  with  results  from  visualization  experiments.  Hartnig 
et  al.  [41  ]  show  a  tomogram  of  a  water  filled  sample  of  a  GDL.  The 
three  different  phases,  the  fibre  material,  the  water  filled  pores  and 
the  free  pores  are  presented  together,  and  additionally  every  phase 
is  shown  in  an  extra  picture.  Regarding  the  water  filled  pores,  one 
can  see  many  small  water  clusters  formed  within  the  GDL  material. 

Regarding  further  experimental  results  [37,38],  one  can  find 
water  cluster  formation  below  the  ribs  within  the  GDL  material. 
However,  one  has  to  take  into  account  that  this  experimental  visu¬ 
alization  of  the  water  within  a  GDL  was  done  for  a  new  GDL 
material  at  different  current  densities.  The  higher  the  current  den¬ 
sity,  the  more  water  clusters  are  formed.  As  within  this  study,  we 
are  focussing  on  variation  of  the  PTFE  content,  we  cannot  directly 
compare  our  results  with  these  experimental  ones.  Nevertheless, 
one  can  conclude  that  our  model  and  the  experiments  agree  insofar 
that  water  cluster  formation  starts  below  the  ribs. 

A  further  aspect  which  can  be  considered  employing  the 
Hoshen-Kopelman  analysis  is  the  distribution  of  cluster  sizes  over 
the  steady-state  simulations.  It  can  be  calculated  from  the  clus¬ 
ter  map  data  which  contains  for  each  cluster  the  respective  size 
averaging  over  different  configurations.  Due  to  the  quite  high  com¬ 
putational  effort  for  the  Hoshen-Kopelman  algorithm,  this  could 
only  be  done  30  times  within  each  steady-state  run  consisting  of 
3  x  106  iterations.  Nevertheless,  the  data  should  be  representative 
enough  to  provide  relevant  information. 

Fig.  7  features  the  corresponding  cluster  size  distribution  den¬ 
sity  curves  for  the  different  PTFE  coverage  values.  There  are 


two  characteristic  trends  which  can  be  monitored  here:  First, 
the  maximum  cluster  size  changes  drastically  with  decreasing 
hydrophobicity.  This  is  especially  evident  between  65%  PTFE  and 
55%  PTFE,  where  the  increase  in  size  amounts  one  order  of  mag¬ 
nitude,  from  103  to  104.  This  corresponds  well  to  the  volatile 
behaviour  between  these  two  ageing  states  as  depicted  above. 


Fig.  7.  Cluster  size  distributions  for  different  PTFE  coverage  values.  The  cluster  size 
is  given  in  voxels. 
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Second,  in  the  size  region  between  approximately  10  and  100 
voxel,  there  is  a  clear  trend  to  increase  the  number  of  these  mid-size 
clusters  going  from  85%  PTFE  to  lower  contents.  The  only  exception 
is  55%  PTFE,  where  the  very  large  cluster  may  ‘absorb’  many  of  the 
smaller  ones. 

4.  Conclusions  and  outlook 

As  the  properties  of  the  water  household  and  the  degradation 
mechanisms  of  fuel  cell  GDLs  are  still  not  well  understood,  there 
is  a  strong  need  also  for  theoretical  models  suitable  to  tackle  this 
problem.  The  MC  model  presented  in  this  work  provides  a  useful 
tool  to  investigate  the  relevant  processes  concerning  the  distribu¬ 
tion  of  liquid  water  within  the  GDL  on  the  p,m  scale.  Although  there 
are  well-established  models  in  the  field,  it  can  serve  as  a  valuable 
completion  due  to  its  special  features  and  the  comparatively  low 
computational  effort. 

The  results  presented  in  this  work  show  a  clear  influence  of 
the  PTFE  coverage  upon  the  amount  of  water  residing  within  the 
GDL  fibre  structure  and  the  structural  properties  of  the  water  phase 
observed.  First,  with  decreasing  PTFE  coverage  more  and  also  larger 
water  clusters  are  formed,  but  upon  reaching  a  critical  value,  the 
behaviour  changes  to  the  formation  of  a  comparatively  very  tall 
cluster  geometrically  connecting  all  GDL  regions.  This  may  be  crit¬ 
ical  for  fuel  cell  operation,  and  could  on  the  other  hand  further 
accelerate  GDL  ageing  in  the  sense  of  a  self-enhancing  mechanism. 

Experimental  validation  of  the  results  can  be  directly  achieved 
by  comparison  to  in  situ  observations  using,  e.g.  synchrotron 
tomography  or  other  imaging  techniques,  but  unfortunately  to  the 
best  of  our  knowledge,  no  appropriate  in  situ  data  on  aged  material 
is  available  for  comparison.  First  appraisal  of  the  results  appears 
promising,  and  in-deep  comparison  is  planned  for  the  future. 
Furthermore,  simulations  employing  GDL  structures  obtained  by 
reconstruction  of  tomography  data  and  larger  simulation  cells  are 
going  to  be  performed  to  reach  further  improvement  in  the  descrip¬ 
tion  of  real  systems. 
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